subroutine rks3
use var
implicit none

    integer i,m,j
    real(8) du,dv
    real(8) Qmean_n_k(1:gridsx,1:gridsy,5)
    
    
    !do KRK=1,3  
    !        call comput_flux1
    !        call comput_flux2
    !        call time_advance           
    !        call comput_original_var 
    !        !if( tstep==1 .and. krk==2) then
    !        !        write(*,*)u(:,1)!,Qmean(i,j,2),Qmean(i,j,4)
    !        !endif
    !        call boundary             
    !    end do
    
    !k=1
    
    call comput_flux1
    call comput_flux2    
    call rhs(Qmean,Qmean_n_k)
    do j=1,gridsy
        do i=1,gridsx
            Qmean(i,j,1)=Qmean_n(i,j,1)-dt*Qmean_n_k(i,j,1)
            
            Qmean(i,j,2)=Qmean_n(i,j,2)-dt*Qmean_n_k(i,j,2)
            
            Qmean(i,j,3)=Qmean_n(i,j,3)-dt*Qmean_n_k(i,j,3)
            
            Qmean(i,j,4)=Qmean_n(i,j,4)-dt*Qmean_n_k(i,j,4)
            !do m=1,4
            !    du=-(flux1(i,j,m)-flux1(i-1,j,m))/dx
            !    dv=-(flux2(i,j,m)-flux2(i,j-1,m))/dy
            !    Qmean(i,j,m)=Qmean_n(i,j,m)+dt*du+dt*dv
            !end do            
            !Qmean(i,j,5)=Qmean_n(i,j,5)-dt*Qmean_n_k(i,j,5)        
        enddo         
    end do
    call comput_original_var
    call boundary
    
    !k=2
    call comput_flux1
    call comput_flux2    
    call rhs(Qmean,Qmean_n_k)
    do j=1,gridsy
        do i=1,gridsx
            Qmean(i,j,1)=3.d0/4.d0*Qmean_n(i,j,1)+1.d0/4.d0*(Qmean(i,j,1)-dt*Qmean_n_k(i,j,1))
            
            Qmean(i,j,2)=3.d0/4.d0*Qmean_n(i,j,2)+1.d0/4.d0*(Qmean(i,j,2)-dt*Qmean_n_k(i,j,2))
            
            Qmean(i,j,3)=3.d0/4.d0*Qmean_n(i,j,3)+1.d0/4.d0*(Qmean(i,j,3)-dt*Qmean_n_k(i,j,3))
            
            Qmean(i,j,4)=3.d0/4.d0*Qmean_n(i,j,4)+1.d0/4.d0*(Qmean(i,j,4)-dt*Qmean_n_k(i,j,4))
            !do m=1,4
            !    du=-(flux1(i,j,m)-flux1(i-1,j,m))/dx
            !    dv=-(flux2(i,j,m)-flux2(i,j-1,m))/dy
            !    Qmean(i,j,m)=3.d0/4.d0*Qmean_n(i,j,m)+1.d0/4.d0*(Qmean(i,j,m)+dt*du+dt*dv)
            !end do
            !Qmean(i,j,5)=3.d0/4.d0*Qmean_n(i,j,5)+1.d0/4.d0*(Qmean(i,j,5)-dt*Qmean_n_k(i,j,5))       
        enddo         
    end do
    call comput_original_var
    call boundary
    
    !k=3
    call comput_flux1
    call comput_flux2
    call rhs(Qmean,Qmean_n_k)
    do j=1,gridsy
        do i=1,gridsx
            Qmean(i,j,1)=1.d0/3.d0*Qmean_n(i,j,1)+2.d0/3.d0*(Qmean(i,j,1)-dt*Qmean_n_k(i,j,1))
            
            Qmean(i,j,2)=1.d0/3.d0*Qmean_n(i,j,2)+2.d0/3.d0*(Qmean(i,j,2)-dt*Qmean_n_k(i,j,2))
            
            Qmean(i,j,3)=1.d0/3.d0*Qmean_n(i,j,3)+2.d0/3.d0*(Qmean(i,j,3)-dt*Qmean_n_k(i,j,3))
            
            Qmean(i,j,4)=1.d0/3.d0*Qmean_n(i,j,4)+2.d0/3.d0*(Qmean(i,j,4)-dt*Qmean_n_k(i,j,4))
            !do m=1,4
            !    du=-(flux1(i,j,m)-flux1(i-1,j,m))/dx
            !    dv=-(flux2(i,j,m)-flux2(i,j-1,m))/dy
            !    Qmean(i,j,m)=1.d0/3.d0*Qmean_n(i,j,m)+2.d0/3.d0*(Qmean(i,j,m)+dt*du+dt*dv)
            !end do
            !Qmean(i,j,5)=1.d0/3.d0*Qmean_n(i,j,5)+2.d0/3.d0*(Qmean(i,j,5)-dt*Qmean_n_k(i,j,5))     
        enddo         
    end do
    call comput_original_var
    call boundary
    
    
    !call comput_flux1
    !call comput_flux2
    !Qmean_n_k(1:gridsx,1:gridsy,1:5)=Qmean(1:gridsx,1:gridsy,1:5)
    !do j=1,gridsy
    !    do i=1,gridsx
    !        du=-(flux1(i,j,1)-flux1(i-1,j,1))/dx
    !        dv=-(flux2(i,j,1)-flux2(i,j-1,1))/dy
    !        Qmean(i,j,1)=Qmean_n(i,j,1)+dt*du+dt*dv
    !        
    !        du=-(flux1(i,j,2)-flux1(i-1,j,2))/dx
    !        dv=-(flux2(i,j,2)-flux2(i,j-1,2))/dy
    !        Qmean(i,j,2)=Qmean_n(i,j,2)+dt*du+dt*dv+dt*Qmean_n_k(i,j,1)*gx
    !        
    !        du=-(flux1(i,j,3)-flux1(i-1,j,3))/dx
    !        dv=-(flux2(i,j,3)-flux2(i,j-1,3))/dy
    !        Qmean(i,j,3)=Qmean_n(i,j,3)+dt*du+dt*dv+dt*Qmean_n_k(i,j,1)*gy
    !        
    !        du=-(flux1(i,j,4)-flux1(i-1,j,4))/dx
    !        dv=-(flux2(i,j,4)-flux2(i,j-1,4))/dy
    !        Qmean(i,j,4)=Qmean_n(i,j,4)+dt*du+dt*dv+dt*(Qmean_n_k(i,j,2)*gx+Qmean_n_k(i,j,3)*gy)
    !        !do m=1,4
    !        !    du=-(flux1(i,j,m)-flux1(i-1,j,m))/dx
    !        !    dv=-(flux2(i,j,m)-flux2(i,j-1,m))/dy
    !        !    Qmean(i,j,m)=Qmean_n(i,j,m)+dt*du+dt*dv
    !        !end do
    !        du=-(gamma1(i,j,1)-gamma1(i-1,j,1))/dx
    !        dv=-(gamma2(i,j,1)-gamma2(i,j-1,1))/dy
    !        Qmean(i,j,5)=Qmean_n(i,j,5)+dt*du+dt*dv            
    !    enddo         
    !end do
    !call comput_original_var
    !call boundary
    !
    !!k=2
    !call comput_flux1
    !call comput_flux2    
    !Qmean_n_k(1:gridsx,1:gridsy,1:5)=Qmean(1:gridsx,1:gridsy,1:5)
    !do j=1,gridsy
    !    do i=1,gridsx
    !        du=-(flux1(i,j,1)-flux1(i-1,j,1))/dx
    !        dv=-(flux2(i,j,1)-flux2(i,j-1,1))/dy
    !        Qmean(i,j,1)=3.d0/4.d0*Qmean_n(i,j,1)+1.d0/4.d0*(Qmean(i,j,1)+dt*du+dt*dv)
    !        
    !        du=-(flux1(i,j,2)-flux1(i-1,j,2))/dx
    !        dv=-(flux2(i,j,2)-flux2(i,j-1,2))/dy
    !        Qmean(i,j,2)=3.d0/4.d0*Qmean_n(i,j,2)+1.d0/4.d0*(Qmean(i,j,2)+dt*du+dt*dv+ dt*Qmean_n_k(i,j,1)*gx)
    !        
    !        du=-(flux1(i,j,3)-flux1(i-1,j,3))/dx
    !        dv=-(flux2(i,j,3)-flux2(i,j-1,3))/dy
    !        Qmean(i,j,3)=3.d0/4.d0*Qmean_n(i,j,3)+1.d0/4.d0*(Qmean(i,j,3)+dt*du+dt*dv+ dt*Qmean_n_k(i,j,1)*gy)
    !        
    !        du=-(flux1(i,j,4)-flux1(i-1,j,4))/dx
    !        dv=-(flux2(i,j,4)-flux2(i,j-1,4))/dy
    !        Qmean(i,j,4)=3.d0/4.d0*Qmean_n(i,j,4)+1.d0/4.d0*(Qmean(i,j,4)+dt*du+dt*dv+ dt*(Qmean_n_k(i,j,2)*gx+Qmean_n_k(i,j,3)*gy))
    !        !do m=1,4
    !        !    du=-(flux1(i,j,m)-flux1(i-1,j,m))/dx
    !        !    dv=-(flux2(i,j,m)-flux2(i,j-1,m))/dy
    !        !    Qmean(i,j,m)=3.d0/4.d0*Qmean_n(i,j,m)+1.d0/4.d0*(Qmean(i,j,m)+dt*du+dt*dv)
    !        !end do
    !        du=-(gamma1(i,j,1)-gamma1(i-1,j,1))/dx
    !        dv=-(gamma2(i,j,1)-gamma2(i,j-1,1))/dy
    !        Qmean(i,j,5)=3.d0/4.d0*Qmean_n(i,j,5)+1.d0/4.d0*(Qmean(i,j,5)+dt*du+dt*dv)       
    !    enddo         
    !end do
    !call comput_original_var
    !call boundary
    !
    !!k=3
    !call comput_flux1
    !call comput_flux2
    !Qmean_n_k(1:gridsx,1:gridsy,1:5)=Qmean(1:gridsx,1:gridsy,1:5)
    !do j=1,gridsy
    !    do i=1,gridsx
    !        du=-(flux1(i,j,1)-flux1(i-1,j,1))/dx
    !        dv=-(flux2(i,j,1)-flux2(i,j-1,1))/dy
    !        Qmean(i,j,1)=1.d0/3.d0*Qmean_n(i,j,1)+2.d0/3.d0*(Qmean(i,j,1)+dt*du+dt*dv)
    !        
    !        du=-(flux1(i,j,2)-flux1(i-1,j,2))/dx
    !        dv=-(flux2(i,j,2)-flux2(i,j-1,2))/dy
    !        Qmean(i,j,2)=1.d0/3.d0*Qmean_n(i,j,2)+2.d0/3.d0*(Qmean(i,j,2)+dt*du+dt*dv+ dt*Qmean_n_k(i,j,1)*gx)
    !        
    !        du=-(flux1(i,j,3)-flux1(i-1,j,3))/dx
    !        dv=-(flux2(i,j,3)-flux2(i,j-1,3))/dy
    !        Qmean(i,j,3)=1.d0/3.d0*Qmean_n(i,j,3)+2.d0/3.d0*(Qmean(i,j,3)+dt*du+dt*dv+ dt*Qmean_n_k(i,j,1)*gy)
    !        
    !        du=-(flux1(i,j,4)-flux1(i-1,j,4))/dx
    !        dv=-(flux2(i,j,4)-flux2(i,j-1,4))/dy
    !        Qmean(i,j,4)=1.d0/3.d0*Qmean_n(i,j,4)+2.d0/3.d0*(Qmean(i,j,4)+dt*du+dt*dv+ dt*(Qmean_n_k(i,j,2)*gx+Qmean_n_k(i,j,3)*gy))
    !        !do m=1,4
    !        !    du=-(flux1(i,j,m)-flux1(i-1,j,m))/dx
    !        !    dv=-(flux2(i,j,m)-flux2(i,j-1,m))/dy
    !        !    Qmean(i,j,m)=1.d0/3.d0*Qmean_n(i,j,m)+2.d0/3.d0*(Qmean(i,j,m)+dt*du+dt*dv)
    !        !end do
    !        du=-(gamma1(i,j,1)-gamma1(i-1,j,1))/dx
    !        dv=-(gamma2(i,j,1)-gamma2(i,j-1,1))/dy
    !        Qmean(i,j,5)=1.d0/3.d0*Qmean_n(i,j,5)+2.d0/3.d0*(Qmean(i,j,5)+dt*du+dt*dv)     
    !    enddo         
    !end do
    !call comput_original_var
    !call boundary
    
    
    !do j=1,gridsy
    !    do i=1,gridsx
    !        du=-gamma1(i,j,1)/dx
    !        dv=-gamma2(i,j,1)/dy
    !        if(KRK.eq.1)then
    !            Qmean(i,j,5)=Qmean_n(i,j,5)+dt*du+dt*dv
    !        else if(KRK.eq.2)then
    !            Qmean(i,j,5)=3.d0/4.d0*Qmean_n(i,j,5)+1.d0/4.d0*(Qmean(i,j,5)+dt*du+dt*dv)
    !        else
    !            Qmean(i,j,5)=1.d0/3.d0*Qmean_n(i,j,5)+2.d0/3.d0*(Qmean(i,j,5)+dt*du+dt*dv)
    !        end if         
    !    enddo
    !     
    !end do
    
    
    
    

end subroutine 